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We present a finite temperature Monte Carlo study of the XY-model in the vortex representation, 
and study its dynamical critical behavior in two limits. The first neglects magnetic field fluctua- 
tions, corresponding to the absence of screening, which should be a good approximation in high 
T c superconductors (k — * oo) except extremely close to the critical point. Here, from finite size 
scaling of the linear resistivity we find the dynamical critical exponent of the vortex motion to be 
z ~ 1.5. The second limit includes magnetic field fluctuations in the strong screening limit [k — » 0) 
corresponding to the true asymptotic inverted XY critical regime, where we find the unexpectedly 
large value z ~ 2.7. We compare these results, obtained from dissipative dynamics in the vortex 
representation, with the universality class of the corresponding model in the phase representation 
with propagating (spin wave) modes. We also discuss the effect of disorder and the relevance of our 
results for experiments. 

PACS numbers: 74.60.-w (Type-Il Sup.), 75.40.Mg (Num. Simulations), 05.70.Fh (Phase Trans.) 



I. INTRODUCTION 

Considerable confusion currently exists, both theoreti- 
cally and experimentally, regarding the dynamical uni- 
versality class of the zero-field superconducting phase 
transition in high temperature superconductors. The 
short coherence length in these materials leads to large, 
non-gaussian fluctuations, and there is some experimen- 
tal evidence that the s±atic behavior is that of the 3- 
dimensional XY modelErO, as in the lambda transition 
in superfluid helium. The dynamical universality class of 
the lambda transition in superfluid helium is that of a 
two component order parameter coupled to a conserved 
density which gives a propagating mode (second sound 
in 4 He; a spin wave in the XY spin model) in the broken 
symmetry phase. Thisjs model E in the notation of Ho- 
henberg and Halperinjj However, the dynamical behav- 
ior of a superconductor could be very different, since the 
lattice acts as a momentum sink destroying Galilean in- 
variance making the system more like helium in a porous 
medium. The existence of disorder does not destroy the 
Goldstone mode in the ordered phase but it does affect 
the vortices and the normal-fluid quasiparticles which 
tend to equilibrate with the lattice rather than co-moving 
with the condensate. 

Furthermore, in d-wave superconductors, scalar dis- 
order causes pair-breaking and quasiparticle branch re- 
combination which may make it inappropriate to assume 
particle number conservation in the dynamics. There 
does not appear to have been any work to date on this 
question. 

A related issue is that of the role of the long-range 
Coulomb interaction, which severely suppresses longitu- 
dinal current fluctuations, leaving only the transverse 
currents associated with vortices. However it is also pos- 
sible that the low superfluid density and screening from 



the high normal fermion density will turn on thc-Jow en- 
ergy longitudinal Carlson-Goldman fluctuationsE3 of the 
order parameter as T c is approached. These microscopic 
fermionic effects may further confuse the data analysis 
and affect the width of the critical regime. 

The issues raised above remain largely unresolved. The 
particular issues that we will explicitly discuss here are 
the role of magnetic screening and the role of disorder. 
In an extreme type-II, system, coupling to gauge fluctu- 
ations is weakjd but nonetheless, in principle, magnetic 
screening becomes important extremely close to the crit- 
ical point-where the system crosses over to inverted XY 
behavior lTjii3 The static correlations related to this is- 
sue haste recently been studied numerically by Olsson and 
TeitelEJ Here we will address the dynamics. 

Most theoretical studies of critical dynamics in XY like 
spin systems have focussed on Landau-Ginsburg repre- 
sentations of the problem involving the phase (i.e. angle) 
of the spin. For static pro pert ies, there also exist equiv- 
alent dual representationsE3li§iLZl in terms of interacting 
vortex degrees of freedom. Although the static proper- 
ties of the phase and vortex representations are the same, 
there is no reason, a priori, why the dynamical univer- 
sality classes should be the same. 

For the reasons discussed above it may be more ap- 
propriate, for superconductors, to consider a model with 
overdamped dissipative dynamics of the topological de- 
fects (the vortices). This is the approach that we will 
take here. 

To add dynamics to the the phase representation, one 
can either include just dissipative dynamics, model A,EI 
in which case the dynamical exponent, z, is close to 2, or 
one can incorporate the propagating (spin wave) modes, 
model E,EI for which z is exactly 3/2 in three dimen- 
sions, (d/2 in d-dimensions) . For the vortex representa- 
tion, which has discrete variables, the natural dynamics 
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is purely dissipative, such as that generated by Monte 
Carlo simulations (in which Monte Carlo time is equated 
with real time). Naively, it would seem unlikely that the 
dissipative dynamics of the vortex representation would 
be in the same universality class as the dynamics of model 
E (phase representation), which has propagatina-modes. 
Surprisingly, recent results by Weber and Jensenlia come 
to the opposite conclusion. They find, for unscreened 
vortex interactions (k = oo) and overdamped dynamics, 
that the dynamical exponent is z — d/2 — 1.5, precisely 
the value one expects in model E dynamics, and consider- 
ably less than the value generally found with dissipative 
dynamics (z w 2). 

In this paper we present results of Monte Carlo cal- 
culations of the dynamical critical exponents for the 
3-dimensional XY model, in the vortex representation, 
with and without magnetic screening. For no screen- 
ing, we., confirm the unexpected result of Weber and 
JensenEfl, while by contrast, for strong magnetic screen- 
ing, we find a rather large enhancement of z. Al- 
though it is-kpown from the Harris criterion and verified 
numericallyo that uncorrelated disorder is weakly irrele- 
vant at the 3-dimensional XY critical point, the effect p£ 
such disorder on the dynamical properties is unknown.EH 
We therefore also investigate the effects of disorder on 
the model with strong screening. 



II. THE MODELS 

The model under investigation here is the XY-model 
with a fluctuating vector potential, 
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where J is the coupling constant (set to unity in the 
simulations) , the <pi denote the phases of the condensate 
on the sites i of a simple cubic lattice of size N = L 3 
with periodic boundary conditions. The sum is taken 
over all nearest neighbors Additionally, we have a 

fluctuating gauge potential ay with the gauge constraint 
[V-a] = 0, and Ao denotes the bare screening length. The 
last term describes the magnetic energy, where the sum 
runs over all elementary plaquettes on the lattice and the 
curl is the directed sum of the gauge potential round a 
plaquette. 

In our work, the main focus is equilibrium vortex dy- 
namics. However, vortices are only implicitly present in 
the above model through the relation 



<j> V0(r) • dr = 27m, 



(2) 



where n = 0,±1,±2,... denotes the net vorticity encir- 
cled by the integration contour. In order to analyze the 
critical behavior of this model due to vortex fluctuations 



explicitly, it is easier to go from the above phase repre- 
sentation of the XY-model to its vortex representation. 
This is achieved by replacing the cosine in Eq. (|l|) with 
the periodic Villain, fun ction and performing fairly stan- 
dard manipulationsliSlijiiZl to obtain 



rii ■ rij Gij[X ]. 



(3) 



Here, the are vortex variables which sit on the links of 
the dual lattice (which is also a simple cubic lattice) and 
Gij is the screened lattice Green's function 



GV/[A ] — J 



L 3 



exp[ik ■ (rj - rj)] 
2 EL=iI 1 -cos(k m )} + A 



(4) 



In the long range case, Ao — > oo, the divergent k = 
contribution has to be excluded from the sum and the 
constraint Yli n i = nas to be imposed. In the short 
range case (finite Ao) this is not necessary. The transfor- 
mations involved in going from the phase to the vortex 
representation also yield the local constraints [V-n]j = 0, 
i.e., there are no magnetic monopoles (zero divergence 
constraint). It is important to note, that the vortex rep- 
resentation does not contain spin wave degrees of free- 
dom anymore, since they have been "integrated out" in 
the transformation procedure. 

We will be interested in two limits of the above model: 
(i) no screening (Ao — > oo), corresponding to the ex- 
treme type-II limit (k — > oo), where the individual vortex 
lines have long range interactions; (ii) strong screening 
(Ao — > 0), i.e., short range interactions, which is sup- 
posed to be the correct description of a [Superconduc- 
tor extremely close to the critical point, El though the 
size of the critical region where such screening is rele- 
vant may be too small to be observable in practice. In 
this limit the interaction reduces to Ga — J(27rAo) 2 and 
Gi^j = (plus exponentially small corrections of or- 
der exp(— r/Ao)). We will in this case use units where 
J(27rAo) 2 = 1. The resulting Hamiltonian is of the very 
simple form 



(Ao 0). 



(5) 



Note, however, that this Hamiltonian is not trivial, since 
the constraints on the local divergence effectively gen- 
erate interactions among the n^. Note further, that 
this is also the dual representation of the XY model 
without a cre eping, in which the temperature scale is 
invertedOO'Eil The static universality class of Eq. (||) 
is then the same as that o£ Eq. (|J) with Ao = oo, and 
is given by XY exponents. B3 The dynamical universality 
class, however, may be different and determining it is one 
of the goals of our study. 

In addition to the pure short-range model, we also 
study the short-range model with quenched random local 
T c . In this case we replace Eq. (||) by 
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Hv = ^X)^ n ^' ( A ° ~* °> &p random )- (6) 



where ^ is uniformly distributed in the interval 
[0.5,1.51. 



III. MONTE CARLO SIMULATION AND FINITE 
SIZE SCALING 

We simulate the following model Hamiltonians in the 
vortex representation: 



1. Eq. 



with An 



2. Eqs. (0) and (0), which corresponds to An 
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We take simple cubic lattices of size L 3 where 4 < L <12- 
64. Periodic boundary conditions are imposed. We start 
with configurations with all = 0, which clearly satis- 
fies the constraints, and a Monte Carlo (MC) move con- 
sists of tradng to create a closed vortex loop around a 
plaquetteBj This trial state is accepted according to the 
heat bath algorithm with probability 1/[1 + exp(/3A2S)], 
where AE is the change of energy and f3 = 1/T. 

Each time a loop is formed it generates a voltage pulse 
AQ = ±1 perpendicular to its plane, the sign depending 
on the orientation of the loop. This leads to a net electric 
fieldEil 



E(t) = ^J v (t) 



with J V (t) = ^, (7) 



where J v (t) is the vortex current density, and At = 1 for 
one full sweep through the system, where, on average, an 
attempt is made to create or destroy one vortex loop per 
plaquette. 

The nonlinear I-V characteristics of the inverted XY 
model can be modeled as the electric field E, due to vor- 
tex current response in the presence of a uniform Lorentz 
force on the vortex lines, proportional to the applied cur- 
rent density JEB In addition, the linear response resis- 
tance can be calculated from the equilibrium voltage- 
voltage fluctuations via the Kubo formulac3 
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R = ^f E A *<n*)n°)> 
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Here, (• • ■) denotes the thermal average, and the voltage 
across the sample is V(t) = LE(t). The resistivity is 
p = L d - 2 R. 

Since we are working with lattices of finite length L, 
one has to employ finite size scaling techniques to extract 
the critical behavior. A detailed scaling theory-has been 
developed for superconductors by Fisher et aljj and we 
now summarize the results from it which will be needed 
for our data analysis. 

Near a second order phase transition the linear resis- 
tivity obeys the scaling law 
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where v is the correlation length exponent, z is the dy- 
namical exponent and p is a scaling function. At the 
critical temperature, p(0) becomes a constant and there- 
fore pu n (T c ,L) ~ L~ {2 ~ d+Z \ If we plot the ratio of pn n 
for different system sizes against T, then 



Hp Vm (L)/p Vln (L')} 
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i.e., all curves for different pairs (L,L') should intersect 
and one can read off the values of T c and z. We will re- 
fer to this kind of data plot as the intersection method. 
With the values of T c and z determined by the intersec- 
tion method we can then use a scaling plot according to 
Eq. (Q) to obtain the .value of v. 

A similar analysisErE3 shows that, above a characteris- 
tic current scale J n \, which varies as ~ L~ 2 at the critical 
point, non- linear response sets in and the electric field 
varies as 
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It is useful to locate the critical temperature from equi- 
librium properties instead of the dynamic scaling of the 
linear resistivity since such simulations are easier to con- 
verge. In the short-range case we do this by considering 
an ensemble with a fluctuating winding number W , de- 
fined by 



W * = L 
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From the finite size scaling relation 

(W 2 ) = /(LV"(T - T c )) 



(12) 
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we can locate the critical point in the short range case, 
both with and without disorder using the fact that the 
winding number is scale invariant at the critical point. 
Note that the simulations of this quantity require global 
moves, where vortex lines going all the way through the 
system are created and destroyed, and therefore do not 
represent the dynamics of the system in a realistic way. 
This does not pose a problem since we use this calculation 
solely as a way to accurately locate the critical point and 
not to follow the dynamics. 

Because different winding number classes are difficult 
to equilibrate efficiently in large systems in this model, 
we use an (exact) duality transformation from the short- 
range vortex models, Eqs. (||) and (^L-jback to an XY 
phase model with a Villain interactionjl3't3 



Zv — ^ <5v-ni,oe 
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J2 e-EiM^^-^- 27 ™ 1 **) (14) 



3 



(a constant prefactor was suppressed in the last equality) . 
Performing the sum over the integer dummy variables 
TTiift leaves a phase-only model. This phase representa- 
tion allows us to take advantage of the Wolff algorithm 
in the simulation which largely circumvents problems of 
critical slowing do-ap and equilibration in the limit of 
large system sizes. E3 The spin- wave stiffness of the dual 
model is given by p s = d 2 f /dA? \^ =0 , where is a con- 
stant vector potential added to the phase gradient in Eq. 
(|l4|). This is related to the winding number fluctuations 
in the vortex model by p s = TL 2 ~ d (W 2 ). 

IV. DATA ANALYSIS 

In this section we analyze our simulation data for the 
XY-model in the vortex representation, starting with the 
model in Eq. (|j|) with long range interactions, i.e., ne- 
glecting screening (Aq — > oo). 

A. Long range interactions 

In Fig. (jl]) we show the data for the linear resistiv- 
ity Pun plotted vs. T obtained from a simulation of the 
Hamiltonian (^|) with Xq = oo. One observes that at 
high temperatures (T = 3.1) there is hardly any size 
dependence in the data, consistent with a correlation 
length which is smaller than the system size. As one 
goes to lower temperature, the size dependence becomes 
stronger. This indicates critical behavior, since in the 
limit L — > oo the linear resistivity should go to zero be- 
low T c . 
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FIG. 1. Linear resistivity pii n vs. T for the vortex model of 
Eq. Q without screening, i.e. Ao = oo. 

However, since it is impossible to locate the critical 
temperature by this kind of plot, we show, in Fig. (||), 
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FIG. 2. Linear resistivity, plotted according to the intersec- 
tion method, vs. T with the data of Fig. (Q). At the intersec- 
tion one can read off T c ps 3.01 and y «s —0.45, corresponding 
to 2 w 1.45(±0.05). 

the same data of Fig. ([[]) plotted according to the inter- 
section method. The curves intersect at approximately 
T c = 3.01 (±0.01) and y w -0.45 corresponding to 
z = 1.45 (±0.05), very similar to the— result of Weber 
and Jensen, who find z = 1.51 (±0.03).Ea Also, the value 
of T c agrees nicely with the one obtained earlier from 
simulations by Dasgupta and Halperin.E3 

Having established T c and z we can now perform a 
scaling plot of the data according to Eq. (||). In Fig. (^) 
we plot pi in L 2 ~ d+z vs. L 1 /" (T-T c ) and find that the data 
collapses best with T c = 3.01 ± 0.01, z = 1.5 ± 0.05 and 
v = 0.66 ±0.01. This independent result confirms that 
the dynamical scaling ansatz for py ln yields the expected 
value of the correlation length exponent v as well as a 
consistent value for z. 



B. Short-Range Interactions 

In this section we will first describe how we locate the 
critical point of the pure and disordered short range mod- 
els defined by Eqs. (||) and (jfy, using finite size scaling 
analysis of Monte Carlo data for static quantities. Then, 
using dynamic scaling analysis we determine the dynam- 
ical exponents both from equilibrium vortex dynamics 
simulations, and from driven nonequilibrium simulations. 

In our simulations of Eq. (^), we used 10 6 — 10 7 sweeps 
to calculate averages, and discard the initial 10% of the 
data for equilibration. For the disordered case we average 
over 10 — 100 samples to obtain small fluctuations. 

As noted earlier, for the purpose of locating the crit- 
ical point as precisely as possible for both the pure and 
random short-range models, we found it convenient to 
perform an exact transformation from the inverted XY 
model back to the phase representation. We compute the 
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FIG. 3. Scaling plot of pu n according to Eq. ([]). Best scal- 
ing was achieved with T c = 3.01 (±0.01), v = 0.66 (±0.01), 
and z = 1.5 (±0.05). 



spin stiffness p s using the Wolff algorithm to overcome 
the critical slowing down and the difficult problem of 
equilibrating different winding number classes.EZEj From 
hyperscaling, jAF 2 } = L d ~ 2 p s /T is scale invariant at the 
critical point £1 Thus curves for different L all cross at 
T = T c as shown in Fig. ^. Using this method, we find for 
the pure model T c = 0.333 ± 0.001 (which agrees nicely 
with the inverse of T c for the long range model) and 
T c = 0.313 ± 0.001 for the disordered model. Of course, 
this technique cannot be used to speed up the determi- 
nation of z since the Wolff algorithm intentionally does 
not represent local relaxation dynamics. Note also that 
the values of (W 2 ) at the crossing points in Fig. ^ in 
the pure (a) and disordered (b) cases agree within the 
error bars, whifib|is what we expect from two-scale fac- 
tor universality^ and the fact that disorder is weakly 
irrelevant to the statics. 

Simulations were then carried out at the measured crit- 
ical temperatures to compute the dynamics in the vortex 
representations. In Fig. (|) the dynamical exponent z of 
the 3D loop model is determined from Monte Carlo data 
for the linear ptesistivity computed from the Kubo for- 
mula, Eq. d)£3 Using p lin - L 1 - 2 at T = T c , and from 
a power law fit to the data we get z — 2.7±0.1. This value 
is consistent with, but more accurate than, the result of 
Wengel and Youngcil. We also verified that is possible 
to collapse the data for different sizes and temperatures 
using Eq. (g) with v w 2/3. 

In the case of the disordered vortex model, it is con- 
venient to determine z from the nonlinear I-V character- 
istics instead of the linear resistivity. Driving the sys- 
tem out of equilibrium makes it cheaper to converge the 
simulations, which is desirable for the disorder averag- 
ing. Fig. (g) shows data for the nonlinear I-V charac- 
teristics. For the largest size studied, L = 64, it is rea- 
sonable to assume that the data is in the range where 
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FIG. 4. Plot of Monte Carlo data for the the winding num- 
ber fluctuations, (W 2 ) vs. temperature T for (a) pure and (b) 
dirty 3D screened vortex models. According to scaling theory 
(see text) the critical temperature is where curves (formed by 
straight-line interpolation between data points) intersect. 



J > J n i ~ l/L 2 . Hence, according to Eq. (pi]), the dy- 
namical exponents can be obtained from power law fits 
of the form E ~ j( 1+2 )/ 2 to the data at T c . The fits were 
done in the current interval where a power law best de- 
scribes the data, leaving out the highest currents where 
saturation artifacts in the simulation limit the voltage. 
This gives z pure = 2.60 ± 0.1 and Zdirty = 2.69 ± 0.1. 
The small discrepancy between these values is within the 
statistical uncertainty of the simulations, and the values 
are consistent with the exponent obtained from the cal- 
culation of the linear resistivity The coincidence of z for 
the pure model with the result of linear resistivity gives 
a consistency test showing that the nonlinear I-V charac- 
teristics correctly determines z. In particular, our model 
assumption of a uniform current driving the vortices in 
the presence of the disorder does not seem to introduce 
any errors. 



V. SUMMARY AND CONCLUSION 

We have performed simulations of the dynamics of 
the 3-dimensional XY model in a vortex representation 
with and without magnetic screening. Without screen- 
ing, we find z m 3/2, in agreement with earlier work of 
Weber and Jensenta, who note that this result agrees 
with the dynamical critical exponent of the phase model 
with spin wave degrees of freedom, model Eli. How- 
ever, the spin wave degrees of freedom are separated 
from the vortex degrees-of freedom when going to the 
vortex representatiorE3~£ZI and the remaining vortex de- 
grees of freedom have only dissipative motion. Hence 
we find it quite surprising, that the exponents from 
spin wave and vortex dynamics agree numerically. We 



5 




FIG. 5. Determination of the dynamical critical exponent 
from MC data for the linear resistivity p computed using the 
Kubo formula for the pure screened vortex model. A power 
law fit to the data (straight line) for T = T c gives the exponent 
z = 2.7 ± 0.1. Inset: Data plotted using the intersection 
method similarly as in Fig. [| The curves have L/L' =4/6 
(triangles up), 6/8 (triangles left), 8/12 (triangles down), and 
12/16 (triangles right) and give the same values for z. 

note that there is some experimental evidence for super- 
diffusive 3-dimensional XY dynamics in superconductors 
with z ~ 1 — 1.5, but only for the case of finite mag- 
netic fields Mil 

One can try to argue that it is the long-range forces 
among the vortices that accounts for the super-diffusive 
(i.e. z < 2) behavior. Indeed one can even argue that 
the spin waves are implicitly present since they are what 
mediate the long-range forces. However, we note that, 
in this model, the long-range forces are instantaneous 
and not retarded. Hence it is still a mystery to us why 
the vortex model appears to be consistent numerical]* 
with Model E dynamics. Note that Lee and Stroud,c£l 
who measured I — V characteristics on the resistively 
shunted junction model (which is described by Model E 
dynamics) using Langevin (as opposed to Monte Carlo) 
dynamics, find z — 1.5 (±0.5), as expected for this model. 

The value of z for the short-range case does not seem 
to be significantly affected by the disorder, and so dis- 
order appears to be irrelevant dynamically as well as 
statically. In both the pure and dirty case however, z 
is significantly larger than is usually seen in relaxational 
dynamics where z is typically only slightly larger than 
twoB There is some experimental evidence for such en- 
hanced values of z. Anlage's groupa finds v = 1.0 ± 0.2, 
and z = 2.65 ±0.3. The value of z is measured directly at 
the critical point from the resistivity scaling and is prob- 
ably more reliable than the value of v which is measured 
somewhat more indirectly (necessarily) using data away 
from the critical point £11 It is possible therefore that the 
value of v is actually consistent with the 3DXY value of 
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FIG. 6. Determination of the dynamical critical exponent 
from MC data for the nonlinear I-V characteristics for system 
size 64 3 . Power law fits (solid lines) in the interval of currents 
where the data is best described by a power law, with slope 
(1 + z)/2 determines z. The dynamical exponents in both 
cases agree (within the statistical uncertainty) with the value 
z » 2.7 found in Fig. |[ 

0.667. Moloni et all find z = 2.3 ± 0.2. 

These experimental values are consistent with the 
value we obtain here. However, in the extreme type-II 
limit, the inverted XY critical regime is expected to be 
very narrow and difficult to access experimentally. It is 
therefore unclear at this point how significant this agree- 
ment is. Further work is needed to estimate more pre- 
cisely the crossover point to inverted XY behavior in real 
materials. 

To conclude, we have raised a number of issues 
about experimental and theoretical uncertainties regard- 
ing the dynamical universality class of the superconduct- 
ing phase transition in high T c superconductors. Clearly 
considerably more work needs to be done to address this 
problem. One aspect that we have not yet addressed 
is the effect of disorder in the case of long-ranged un- 
screened interactions. A second problem worth pursuing 
is the precise theoretical relationship between the various 
possible spin wave dynamics in the phase representation 
and the corresponding dynamics of the same model in 
the various dual (vortex) representations. 
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3D pure loop model 

from fit: slope=p=1-z=-1 .68 => z=2.7 +/- 0.1 
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